//
// Ephi - simulation of magnetic fields and particles
// Copyright (C) 2007 Indrek Mandre <indrek(at)mare.ee>
// For more information please see http://www.mare.ee/indrek/ephi/
//
// This program is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation; either version 2 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License along
// with this program; if not, write to the Free Software Foundation, Inc.,
// 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
//

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <assert.h>

#include "ephi.hpp"

int main (int argc, char *argv[])
{
  Statics statics;
  make_polywell_cube (statics, 0.15, 0.035, 0.01, 1000000, 0);

  for ( size_t i = 0; i < 2; i++ )
    {
      const transf3d_pair *pairs;
      size_t N;
      if ( i == 0 )
        {
          N = 24;
          pairs = get_cube_polywell_transforms24 ();
        }
      else
        {
          N = 48;
          pairs = get_cube_polywell_transforms48 ();
        }
      char path[128];
      sprintf (path, "test17_%02d.pov", (unsigned int)N);
      FILE *fp = fopen (path, "w");
      assert (fp);
      fprintf (fp, "#include \"scene.inc\"\n");
      fprintf (fp, "camera { location <2, 1.5, 3>*0.25 look_at <0, 0, 0> sky <0, 1, 0> right -x*image_width/image_height }\n");
      fprintf (fp, "\n");

      fprintf (fp, "#declare mypoly = union {\n");
      std::string poly = statics.get_povray_script();
      fprintf (fp, "%s", poly.c_str());
      vect3d pos(0.3, 0.075, 0);
      vect3d vel(0, 0, 1);
      for ( size_t i = 0; i < N; i++ )
        {
          vect3d v = pairs[i].f * vel;
          vect3d p1 = pairs[i].f * pos;
          vect3d p2 = p1 + v.normal() * 0.05;
          fprintf (fp, " cone { <%f, %f, %f>, 0.02, <%f, %f, %f>, 0 texture { ArrowTexture } }\n",
              prec2double (p1.x), prec2double (p1.y), prec2double (p1.z),
              prec2double (p2.x), prec2double (p2.y), prec2double (p2.z)
          );
        }

      pos = vect3d (0.2, 0.2, 0.25);
      vel = vect3d (-1, 1, 0);
      for ( size_t i = 0; i < N; i++ )
        {
          vect3d v = pairs[i].f * vel;
          vect3d p1 = pairs[i].f * pos;
          vect3d p2 = p1 + v.normal() * 0.05;
          fprintf (fp, " cone { <%f, %f, %f>, 0.02, <%f, %f, %f>, 0 texture { ArrowTexture2 } }\n",
              prec2double (p1.x), prec2double (p1.y), prec2double (p1.z),
              prec2double (p2.x), prec2double (p2.y), prec2double (p2.z)
          );
        }

      fprintf (fp, "}\n");

      fprintf (fp, "object { mypoly rotate <0, clock*360, 0> }\n");

      fclose (fp);
    }

  for ( size_t i = 0; i < 30; i++ )
    {
      const transf3d_pair *pairs;
      pairs = get_cube_polywell_transforms48 ();
      char path[128];
      sprintf (path, "test17_48b%02d.pov", (unsigned int)i);
      FILE *fp = fopen (path, "w");
      assert (fp);
      fprintf (fp, "#include \"scene.inc\"\n");
      fprintf (fp, "camera { location <2, 1.5, 3>*0.25 look_at <0, 0, 0> sky <0, 1, 0> right -x*image_width/image_height }\n");
      fprintf (fp, "\n");

      fprintf (fp, "#declare mypoly = union {\n");
      std::string poly = statics.get_povray_script();
      fprintf (fp, "%s", poly.c_str());

      prec_t angle = 2 * PREC_PI / 4;
      prec_t phi = i * angle / 30;

      vect3d pos(0.3, cos(phi) * 0.2, sin(phi) * 0.2);
      vect3d vel(0, -sin(phi) * 0.2, cos(phi) * 0.2);
      for ( size_t j = 0; j < 48; j++ )
        {
          vect3d v = pairs[j].f * vel;
          if ( j & 1 ) v = -v;
          vect3d p1 = pairs[j].f * pos;
          vect3d p2 = p1 + v.normal() * 0.05;
          fprintf (fp, " cone { <%f, %f, %f>, 0.02, <%f, %f, %f>, 0 texture { ArrowTexture } }\n",
              prec2double (p1.x), prec2double (p1.y), prec2double (p1.z),
              prec2double (p2.x), prec2double (p2.y), prec2double (p2.z)
          );
        }

      fprintf (fp, "}\n");

      prec_t fangle = i * 360.0 / 60.0;
      fangle = 0;
      fprintf (fp, "object { mypoly rotate <0, -%f, 0> }\n", fangle);

      fclose (fp);
    }

  std::vector<vect3d> planes24;
  planes24.push_back (vect3d(0, 1, 0));
  planes24.push_back (vect3d(0, 0, 1));
  planes24.push_back (vect3d (1, -1, 0).normal());
  planes24.push_back (vect3d (1, 0, -1).normal());

  for ( size_t i = 0; i < 100000; i++ )
    {
      vect3d pos(prec_t_drand() * 0.3 - 0.15, prec_t_drand() * 0.3 - 0.15, prec_t_drand() * 0.3 - 0.15);
      vect3d vel(prec_t_drand() * 0.3 - 0.15, prec_t_drand() * 0.3 - 0.15, prec_t_drand() * 0.3 - 0.15);
      vel.normalize();
      for ( size_t j = 0; j < 2; j++ )
        {
          vect3d obf, oef, tbf, tef, cbf, cef, rbf, ref, ov, cv, of, cf;

          vect3d p1 = pos;
          ov = vel;

          statics.getFields (p1, obf, oef);
          of = vel % obf;

          const transf3d_pair *t;
          if ( j == 0 )
              t = &cube_polywell_transform24 (p1);
          else
              t = &cube_polywell_transform48 (p1);
          statics.getFields (p1, tbf, tef);
          cf = t->f * of;
          cv = t->f * ov;
          cbf = t->f * obf;
          cef = t->f * oef;
          rbf = t->b * tbf;
          ref = t->b * tef;

#if 0
          pos.print("pos");
          p1.print("p1");
          (t->f * pos).print ("t*pos");
          obf.print("obf");
          rbf.print("rbf");
          tbf.print("tbf");
          cbf.print("cbf");
#endif
          if ( (tbf - cbf).length() > 0.000001 * obf.length() )
            {
              char b1[128];
              fprintf (stderr, "NOK: bfields dont match(1) at position %s\n", pos.sprint(b1));
              break;
            }

          if ( (obf - rbf).length() > 0.000001 * obf.length() )
            {
              char b1[128];
              fprintf (stderr, "NOK: bfields dont match(2) at position %s\n", pos.sprint(b1));
              break;
            }

          vect3d tf = cv % tbf;
#if 0
          of.print ("of");
          tf.print ("tf");
          cf.print ("cf");
#endif

          if ( j == 0 )
            {
              if ( (tf - cf).length() > 0.000001 * of.length() )
                {
                  char b1[128];
                  fprintf (stderr, "NOK: forces dont match(3) at position %s\n", pos.sprint(b1));
                  break;
                }

              for ( size_t k = 0; k < planes24.size(); k++ )
                {
                  if ( p1 * planes24[k] < -1e-23 )
                    {
                      char b1[128];
                      fprintf (stderr, "NOK: translation not confined to sector at position %s\n", pos.sprint(b1));
                      break;
                    }
                }
            }

          const transf3d_pair *pairs;
          size_t N;
          if ( j == 0 )
            {
              N = 24;
              pairs = get_cube_polywell_transforms24 ();
            }
          else
            {
              N = 48;
              pairs = get_cube_polywell_transforms48 ();
            }

          for ( size_t k = 0; k < N; k++ )
            {
              vect3d tp = pairs[k].f * pos;
              vect3d tv = pairs[k].f * vel;
            }
        }
    }

  return 0;
}

